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Abstract 

Two-step hybrid methods specially adapted to the numerical integration of 
perturbed oscillators are obtained. The formulation of the methods is based 
on a refinement of classical Taylor expansions due to Scheifele [Z. Angew. Math. 
Phys., 22, 186-210 (1971)]. The key property is that those algorithms are able to 
integrate exactly harmonic oscillators with frequency uj and that, for perturbed 
oscillators, the local error contains the (small) perturbation parameter as a factor. 
The methods depend on a parameter v = Loh, where h is the stepsize. Based on 
the B2-series theory of Coleman [IMA J. Numer. Anal, 23, 197-220 (2003)] we 
derive the order conditions of this new type of methods. The linear stability and 
phase properties are examined. The theory is illustrated with some fourth- and 
fifth-order explicit schemes. Numerical results carried out on an assortment of 
test problems (such as the integration of the orbital motion of earth satellites) 
show the relevance of the theory. 
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1 Introduction 



In the last decades, there has been a great interest in the research of methods for 
the numerical integration of initial value problems (IVP) associated to second-order 
ordinary differential equations (ODE) 



in which the first derivative does not appear explicitly. These problems appear often 
in practice. Of course, since (11.11) can be written as an IVP for a system of two equa- 
tions of first-order, the problem can be solved by algorithms for first-order equations. 
However, this will be less efficient than if methods specially devised for the given 
problem would be used. The construction of methods specialized for (II. ip is a well 
established area of investigation. Many multistep methods (such as Stormer-Cowell 
methods) and two-step hybrid methods for (II. ip have been developed, see for example 
Lambert & Watson (1976), Chawla (1984), Chawla & Rao (1987), Coleman (1989), 
Simos (1999), Tsitouras (2003), Coleman (2003) and Franco (2006a) to mention a few. 
Two-step hybrid methods are considered to be more efficient than the rival Runge- 
Kutta-Nystrom methods for (II. ip . For example, the standard fourth-order explicit 
Runge-Kutta-Nystrom method (see Hairer et al. (1993)) requires three function eval- 
uations whereas the fourth-order explicit Numerov method of Chawla (1984) requires 
only two function evaluations per step. 

Quite often the solution of (11.10 exhibits an oscillatory behaviour; think, for 
instance, of the pendulum problem in celestial mechanics or of the Schrodinger equa- 
tion in quantum mechanics. For problems having highly oscillatory solutions standard 
methods with unspecialized use can require a huge number of steps to track the oscilla- 
tions. One way to obtain a more efficient integration process is to construct numerical 
methods with an increased algebraic order. On the other hand, the construction and 
implementation of high algebraic order methods is not evident. Alternatively, one can 
consider methods that use the detailed information of the high-frequency oscillation. 
There is a vast literature on this subject; an extensive bibliography is summarized by 
Petzold et al. (1997). Scheifele (1971) was concerned with the solution of perturbed 
oscillators, i.e., second-order problems of the form 



y" = f{x,y) 



y{xo) = yo 



y'{xo) = y'o, 



(1.1) 



y" = -u'^y + g{x,y), 



y{xo) = yo, 



y'ixo) = ?/o. 



(1.2) 
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where the magnitude of the perturbation force satisfies \g{x,y)\ « uP' \y\. Scheifele 
rewrote the solution of f 1 1.2 1) as a series of a set of functions, the G-functions, more 
adequate to perturbed oscillators than the classical polynomial Taylor expansion. The 
Scheifele G-functions method is capable to integrate exactly the harmonic oscillator 
or unperturbed problem (i.e. (11. 2p with (7 = 0). In spite of its excellent behaviour, 
the Scheifele G-functions method has the disadvantage that it is strictly application- 
dependent. Several authors have applied Scheifele's approach for constructing numer- 
ical methods adapted to perturbed oscillators. Most of these papers are focused on 
space dynamical problems such as an accurate integration of orbit problems or long- 
term prediction of satellite orbits. Some Scheifele G-functions based multistep codes 
are designed by Martin & Ferrandiz (1997). Also adapted methods without first deriva- 
tives have been constructed by Lopez et al. (1999). A first Runge-Kutta type version 
of the Scheifele G-functions method is due to Gonzalez et al. (1999). A theoreti- 
cal foundation for these adapted Runge-Kutta-Nystrom (ARKN) methods is given by 
Franco (2002,2005,2006b). 

Our objective in this paper is apply Scheifele's approach to two-step hybrid meth- 
ods. This was already proposed by Van de Vyver (2007a) for the simple explicit Nu- 
merov method. The excellent numerical results reported in that paper strongly suggest 
to construct higher-order methods of this type. This is possible when a more theoret- 
ical framework would be developed. This is the purpose of this work. The paper is 
organized as follows. Section [2] is of an introductory nature: we recall the class of 
classical two-step hybrid (TSH) methods. In Section [3] we recall Scheifele's approach. 
This idea will be extended to TSH methods, the resulting methods are denoted by 
ATSH methods. Section H] is devoted to the order conditions for ATSH methods. This 
part heavily relies on the work of Coleman (2003) for classical TSH methods. Some 
general stability results for ATSH methods are reported in Section [51 The concepts 
of such a stability analysis find its origin in the work of Coleman & Ixaru (1996) and 
Franco (2005). Section [H] provides general results on the phase properties of ATSH 
methods. The analysis is based on the work of Franco (2005). Section [7] deals with the 
construction of fourth- and fifth-order explicit ATSH methods. Several possibilities are 
explored such as minimizing the error constant, increasing the phase-lag order, dissipa- 
tive or not, . . . The classical companions of the new methods are previously derived by 
Franco (2006a). Section [8] collects numerical examples for a variety of problems chosen 
to illustrate particular features of the ATSH methods obtained. The new methods are 
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compared with other high-quahty methods. The paper concludes with a brief summary 
of the work considered here. 



2 Classical two-step hybrid methods 

Two-step hybrid (TSH) methods for (11.11) are defined by 

s 

Yi = (1 + Ci) yn - Ci Vn-i + h^^ aij f{xn + Cj h,Yj), i = 1, . . . , s, (2.3) 



n+l = 2yn- Vn-l + f {Xn + Ci h, Yi] 



(2.4) 



i=l 



where i/n-i, Un and yn+i are approximations of ?/(x„ — h), y{xn) and ?/(x„ + h), re- 
spectively. TSH methods can be in short-hand notation represented by the Butcher 
table 
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where c, 6 G Mf^^ and A G W^^. These coefficients are derived by imposing the 
necessary and sufficient conditions for convergence, i.e. consistency and zero-stability, 
see Henrici (1962) for the general theory. 

For exact starting values, the local truncation error (Ite) of the method at Xn is 

s 

he = y{xn + h) -2y{xn) + y{xn - h) - ^&i/(a;„ + Cih,Yi). (2.5) 

1=1 

The method is of algebraic order p if lie = 0{h^^'^). The principal local truncation 
error (plte) is the leading term of (12. 5p . For a pth-order method this is of the form 



plte 



(p + 2)! 



a{t) {i + {-ir'-b''m)Fmyn,y'n 



(2.6) 



p(t)=p+2 



where a(t), p{t), \E'"(t), F{t) and T2 are defined in Coleman (2003). The coefficients of 
F{t){iin,y'n) i'^ (12. 6p will be denoted as ep+i(t). The quantity 



^2 

p+2 
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\p(t) = 



(2.7) 



(2.9) 



will be called the error constant of the pth-order method. Traditionally, the order con- 
ditions for TSH methods are usually derived by expansions in Taylor series. These ex- 
pansions are calculated essentially by brute force. On the other hand, Coleman (2003) 
obtained the order conditions for TSH methods by using the theory of B-series. Analo- 
gously to the case of RK(N) methods, the determination of the order of a TSH method 
is based on checking certain relationships between the coefficients of the method. 

The linear stability analysis of methods for solving (11.11) is based on the scalar 
test equation (see Lambert & Watson (1976)) 

y" = -X'y, A>0. (2.8) 

An application of a TSH method to (12. 8p yields 

Y = {e + c)u,-cyr.^i-H^AY, H = Xh, 
yn+i = 2yn-yn~i- H'^b^y, 

where Y = (Yi, . . . , Y^)"^ and e = (1, . . . , 1)^ G W^^. Elimination of the vector Y 
from (12. 9p results in the difference equation 

yn+i-SiH^)yn + PiH^)yn-i = 0, (2.10) 

where 

S{H^) = {I + H^A)-'{e + c), 

(2.11) 

P{H^) = 1 - HH'^ (I + Ay c. 
The solution of the difference equation (12.101) is determined by the characteristic equa- 
tion 

- S{H^)^ + P{H^) = 0. (2.12) 

Of particular interest for periodic motion is the situation where the roots of (12.121) 
lie on the unit circle. For example, in celestial mechanics it is desired that numerical 
orbits do not spiral inwards or outwards. This periodicity condition is equivalent to 

P{H^) = 1 and \S{H^)\<2, yHe{0,H^,,), (2.13) 

and the interval (0, -ffper) is called the interval of periodicity. If the necessary condition 
P(H'^) = 1 to have of a non-empty interval of periodicity is not satisfied, we can ask 
when the numerical solution remains bounded. This stability condition is equivalent 
to 

P{H') < 1 and \SiH')\ < l + P{H^), VH E (0,/flJ, 
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and the interval (0, H^tab) called the interval of absolute stability. 

Another related concept, which is important when solving problems of the form (11.11) 
is the phase-lag of the method. In phase analysis one compares the phases of exp(± i H) 
with the phases of the roots of the characteristic equation (12.121) . Following the ap- 
proach of van der Houwen & Sommeijer (1987) for RKN methods, the quantities 

ct>{H) = H - arccos (^-^^^S^y d{H) = 1 - y^P{lP) , (2.14) 

are the phase-lag (or dispersion) and the dissipation (or amplification error), respec- 
tively. The method is said to have phase-lag order q and dissipation order r if 



(l){H) = H'>+^ + 0{H'i+^), d{H) = Cd H'+^ + 0{H'+^). 

The constants and Cd are called the phase-lag and dissipation constants, respectively. 
Methods with d{H) = are zero-dissipative. 

3 Two-step hybrid methods for perturbed oscilla- 
tors 

3.1 Notations and exact solution 

Although, Scheifele's method is based on G-functions, in this paper we consider the 
related 0-functions which are suggested by Franco (2002) for the derivation of the order 
conditions for ARKN methods. The coefficients of Scheifele's G-functions method are 
dependent on the frequency uj and stepsize h. By using the </)-f unctions, the coefficients 
are dependent on only one variable v = uoh. 

The solution of (11.21) can be expressed as 

sin(z/) 1 

y{xn + h) = y{xn) cos(z/) + hy\xn) ^ " / Q^^i vi^)) sin(u; (x„+i - x)) dx. 

^ ^ Jx„ 

(3.15) 

We carry out the change of variable x = Xn + h z in (13.151) and we denote (p{x) = 
g{x,y{x)). Now the exact solution becomes 

, sin(z/) , 2 /"^ / , siniu (1 — z)) 
y{Xn + h) = y{Xn> cos[u) + hy [Xn) \-h / (p{Xn + hz) dz. (3.16) 
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Suppose that the function (f{x) admits an expansion of the form 



ip{xn + hz) = J2h'v^'\xn)^. (3.17) 



We can write that 



/ 7\ / \ / \ 7 // X sin(z/) , 14.2 (i), \ f sin(i/ (1 — 2;)) , 
y{xn + h)= y{xn) cos z/ + hy'{xn) — — + > /i^+V^'H^^n / — -7 dz. 

(3.18) 

Introducing the following notations 



smu^ 



^ sin(z/ (1 — 2;)) 2;-^ 



0o(z/) = cos(z/), 0i(z/) = , (t)j+2{^) = / q-c?2, j > 0, 

(3.19) 

we arrive to the expression of the exact solution of the perturbed problem (11. 2p in 
terms of 0-functions 

00 

y{xn + h)= y^M^) + hy'M^) + ^^'^M (l>M^). (3.20) 

i=o 

Remark that the analytical solution of the harmonic oscillator is approximated exactly 
by the expansion (13.201) . 

Some interesting properties of the ^-functions are listed in the following theorem. 



Theorem 1 1. YmvchAi') = —, j > 0. 
2. The (p-functions can he expressed as 

E(-l)'^7^)' ^-^0' (3-21) 



'-IV I , ... v^^ 



^^2A^) = ^\^^^K^)- l^y--) (2l)! 

yrf-M-g(-i)'piTT^). ^ao- (3.22) 



3. The Taylor series expansions of the (p-functions are 



00 o u 



)=E(-1)'7JTT-TT. ^-^O- (3-23) 



" (2* + j)! 



4. 0j+i(z/) = / cos{u{l - z)) —dz, j>0. 

Jo J- 

5. We have the following recurrence relation 

0,(i/) + z/V,+2M = ^, J>0. (3.24) 

The 0-functions are related to the Scheifele G-functions by Gj{h) = 4>j{i^), j > 0. For 
further details and proofs about G-f unctions, see Scheifele (1971), Fairen et al. (1994) 
and Martin & Ferrandiz (1997). 

According to Theorem [1] (point 1) it is clear that when the frequency — > 
— > 0) the series fl3.20p will become 

y{xn + h)= y{xn) + hy'{xn) + ^ / ■ ^ o^l ^^''"'^^^"^' ^^'^^^ 

i=o ^^ + ^)- 

which is the classical Taylor expansion of the exact solution. Thus Scheifele's se- 
ries fl3.20p is a refinement of the classical Taylor method. 



3.2 Formulation of the method 

An s-stage TSH method (I2.3p -( !2l4l) can be rewritten in the following alternative form 

s 

s 

yn+1 = 2yn- yn-1 + ^"^ ^ K 

1=1 

We can see that k[ are evaluations of the function / at the points Xn + Q /i, where the 
second argument is an approximation to the solution at this point. Then, we have 

s 

y{xn + Cih)K.{l + a) yn - Ci yn-1 + h'^^ ^ij k'j, i = 1, . . . , s. 

i=i 

For perturbed oscillators, i.e. when f{x,y) = —oj'^y + g{x,y), the internal stages can 
be approximated by 

ki = g{xn + Cih,Yi), i = l,...,s, (3.26) 
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where 

s 

Yi = {1 + d) Un - Ci Un-l + X] "'^^ ^'J + ^j)- 

i=i 

The coefficients represent the weights of the quadrature formulas used in the ap- 
proximation of the internal stages. 

The final stage is determined as follows. We can avoid the calculation of the first 
derivative of the solution of (13.161) by adding this expression with positive and negative 
stepsize to get 

y{xn + h)=2 (Po{i^) y{xn) - y{xn -h) + h^ [ ~ I^D) ^(3,^ + hz)dz. (3.27) 



J-i V 

We shall approximate the exact solution by using the quadrature formula 



/•^ sin(z/(l - \z\)) . ^ , ^ , 
"^-1 i=l 



where the fc- values are given by (I3.26p . 

Altogether, we arrive to the following definition. 

Definition 1 An s-stage adapted two-step hybrid (ATSH) method for the numerical 
integration of the IVP ( fi.^) is given by the scheme 

s 

Yi = {1 + d) yn - Ci yn-i + h'^^ ^ij (^~^^ + g{xn + Cj h, Yj)^, 1, • • • , s, 



s 

yn+1 = 2 0o(l^) Vn - yn-l + h'^ ^kg^Xn + Ci h,Yi), 

i=l 

which can be expressed in Butcher notation by the table of coefficients 



(3.28) 
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Remark that when — ^ 0, ATSH methods reduce to classical TSH methods. 

As said, the convergence of a method is covered by consistency and zero-stability. 
The consistency (i.e. algebraic order is at least 1) follows form Section |H The theorem 
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of Ixaru & Rizea (1987) says that any method apphed to y" = with resuhing difference 
equation 

Vn+l + ai{h) + Vn-l = 0, 

is zero-stable if ai{h) = — 2 + 0{h'^), q > 2. Using Theorem [1] (point 3) it is easy to 
see that ATSH methods are zero-stable. 

4 Order conditions for ATSH methods 

Similarly to the classical case, the principal local truncation error (Ite) of the ATSH 
method (13.281) is given by 

He = y{xn + h) -2 (f)o{i^) y{xn) + y{xn - h) - ^^hi g{xn + Q Yi). 

i=l 

The method is of algebraic order p if Ite = 0{h^^'^). Our next aim is to derive order 
conditions for ATSH methods by adapting the recently developed B2-series theory of 
Coleman (2003). In what follows, the reader is referred to that paper for all the defi- 
nitions and notations. The theory of B2-series is applicable only to one-step methods 
so we have to search for a one-step formulation of ATSH methods. A modification of 
Coleman's proofs at several places will deliver the requested order conditions. 

4.1 Adapted B2-series 

Repeated differentiation of (p with respect to the independent variable x gives 

The difference with the classical theory lies in the fact that every elementary differential 
starts with a Frechet-derivative of g instead of /. The following definition explains how 
each elementary differential can be associated with a rooted tree. 

Definition 2 The function G on T2\{0, r'} is defined by 
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1. G{T){y,y')^g. 

2. Ift= [ti,...,t™]2 then 

G{t){y,y') = g^'^\y) (F(t,)(y,y'), . . . , F(tJ(y,y')) , 
where the function F is recursively defined in Definition 3 of Coleman (2003). 

Analogously to the classical theory, it is obvious that 

a{t)Git){y,y'), (4.29) 

p(t)=j+2 

where a{t) represents the number of distinct monotonic labellings of the vertices of 

t e T2. 

B2-series are defined in Definition 4 of Coleman (2003). Here that definition is 
adopted more pertinent for our methods. 

Definition 3 Let (5 a mapping from T2 to R. The adapted B2-series with coefficient 
function (5 is a formal series of the form 

BiM= E -7w^c^it)mGmy,y')- 

teT2\{0,T'} 

Coleman's fundamental lemma is then reformulated for the adapted case as follows. 

Lemma 1 Let B{(3,y) be a classical B2-series. Then h? g{B{(5,y)) is an adapted B2- 
series, 

h'g{B{(3,y))^B{(3",y), 

with 

(5"{0) = (5"{t') = 0, (5"{t) = 2, 
and for all other t — [ti, . . . , tm\2 £ T2, 

m 

p"{t)^pit){p{t)-i)i[m- 

i=l 

The proof is essentially the same as the original proof. 
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4.2 One-step formulation 

By defining := {i/n+i — (po{^)yn)/h the second equation of (13.281) can be expressed 
as a pair of equations 

Vn = (poij^) Vn-l + h Fn-1, 

Fn = M^)Fn^i-uju4>l{u)yn-i + h{}F ®I)g{Y). 
Now, the one-step formulation takes the form 

Un = M{v)un-i + h^{un^i,h), (4.30) 

with 

M{v) = 



-uj u (t)l{iy) (l)o{iy) 
un={ and $(«„_!, h) = ( ^Y)'giY) 



(4.31) 



and Y is defined implicitly by 

Y = {e + c)®yn-c(^yn^i + h^A®I){-uj^Y + g{Y)) 

= (e + (0o(i^) - 1) c) ® yn-i + /i (e + c) (g) + h'^{A0 I) {-cu^ Y + g{Y)). 

(4.32) 

4.3 Order conditions 

The vector m„ is an approximation for z„ = z{xn, h), where 

/ y{x) ^ 

'(^^h)= y{x + h)-<P,{v)y{x) ■ 
\ h I 

For exact starting values, the Ite of the one-step formulation fl4.30p - (l4.32p is 



(4.33) 



dn = Zn- M(Z/) Zn-1 - h ^{z^-i, h) , (4.34) 

with 



^Zn-uh)= I h I , (4.35) 

{b^®I)g{Y) 

where Y is now defined implicitly by 

y = e ® y{xn-i) + (e + c) ® (y(x„) - y{xn-i)) +h^A^ I) {-u^ Y + g{Y)). 
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Definition 4 The ATSH method liS. 28\) is of algebraic order p when dn = 0{h'^'^^). 

We are now ready to present one of the main results of this paper. 

Theorem 2 The ATSH method ^3. 28\) is of algebraic order p if and only if, for trees 
t G T2 

6^<(t) = (l + (-irW)p(t)!0p(i)H, 
for pit) <p+l but not for some trees of order p + 2. 

Proof. Observing fl4.33l) - fl4.35l) we have that the first component of dn is zero. Each 
component of the vector Y can be expanded as a B2-series 

Y,{xr:) = B(^i,,,y{xn))=Y.^^a{t)i,,{t)F{t){yn,y'^). (4.36) 

The coefficients V'i(^) can be generated recursively by formulas (3.6)-(3.7) of Cole- 
man (2003). We substitute the B2-series (I4.36P into the second component of dn and 
we apply Lemma [H An easy calculation gives 

^ I y{.Xn + h)-2 0o(z/) y{xn) + y{xn - k) - k'^ ^bi g{Yi{xn)) 

\ i=l y 

= h[^ E^'Vi'^-'V2.(^) -J2kB{,p':,y{xn)) . (4.37) 
\ j=i i=i J 

With (I4.29P in mind, the left side of (I4.37P becomes 

00 

2E^'Vi'^-'V2,(^)= Yl h''^'^a{t)M^)Gmyn,y'n). (4.38) 



j = l tST2 

p(t) even 

The right side of (14.371) may be written as 

J2 h yn) = Y.^^ "(^) ^'(^) ^(^)(^/- y'n)- (4-39) 

i=l t&T2 P^^>- 

The theorem follows when comparing (I4.38p and (I4.39p . 

The order conditions up to order six are listed in Table [H 

Remark 1 Reconsidering Section 5 of Coleman (2003) it is obvious that, in order to 
reduce the number of order conditions, the simplifying conditions for ATSH methods 
are the same as for classical TSH methods. 
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l^ijOiCittijC^j = 3 04(Z^) 






- , b-; a^A a^u Cu — — - Sa iv\ -\- a (hai z/l 


^78 
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^79 




Eij,fc^i aij Cjajfe = 2 06('^) 


^7,10 




Ei bi ttij ttjk Cfc = 



Table 1: Order conditions 
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4.4 Error analysis 



From the proof of Theorem [2] it follows that the plte of a pth-order ATSH method is 
given by 

^" ' ■ teT2 

p(t)=p+2 



where b^^^'^ and \E'(°^ represents the 6^- and \E'-values of the corresponding classical TSH 
method. The plte of this classical method for (11.11) reads 

plte^'^ = Yl «(^) (l + (-^r" - i^W(l/n,l/;). (4.40) 

p{t)=p+2 

In order to obtain a connection between plte^^^ and plte^^^^ we need a relationship 
between F(t) and (j(t). This can be easily seen as follows. We consider trees in which 
the root starts with a chain of 3 vertices (including the root) having exactly one son. 
We call such a tree a semi-tall tree. We denote by Tg* the set of semi-tall trees. The 
truncated tree t~ of a semi-tall tree t is obtained by deleting the first two vertices. 
Clearly, the number of semi-tall trees of order p + 2 is equal to the number of trees of 
order p. Using the above terminology, it is easy to see that 

F{t){y,y')+uj^F{t-){y,y') iiteT*, 

F{t){y,y') litiT*. 

We conclude with 

plte^^'^ =plte''"+u'j^^^ J2 «W {l + {-ir'-b^'^^<i^^'\t))F{t-){y^,y'J. 

(4.41) 

For the calculation of the error constant, E^J'f^, we have to consider the coefficients of 
F{t){yn,yn) and the coefficients of cj^ F(t-)(y„, y^) in fOB . Observing flCTD -f lCTl ) 
it is clear that 



G{t){y,y') 



rpATSH 
Fp+1 



' ^ ^ 2 li teT*, 



\p(t)=p+2 / 



with U = < (4.42) 



1 if t ^ T*. 
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5 Linear stability analysis 



Linear stability and phase-lag analysis of ATSH methods is also based on the model 
equation (12.81) . However, this equation has to be rewritten in the following appropriate 
form 

y" = -uj^y-ey, a;2 + e>0, (5.43) 

where uj represents an estimation of the dominant frequency A of (12. 8p . and e = — cu^ 
is the error of that estimation. This modified test equation is prompted by the work 
of Franco (2005) for ARKN methods. At the first sight, one should believe that the 
estimated frequency uj should be equal to dominant frequency A. This is generally a 
satisfying approach but in practical applications it is possible to obtain more accurate 
results for different values of A and u. The cubic oscillator 

y" = -y + ey\ 1/(0) = 1, 2/'(0) = 1, 

provides such an example. Although this is a nonlinear problem, for small e-values 
we may apply linear stability analysis, resulting in A = 1. However, Vigo-Aguiar 
et al. (2004) have proved that more accurate results are obtained when selecting 
UJ = VI -0.75e. 

An ATSH method (Km applied to yields 

Y = {e + c)yn-cyn-{iy^ + z)AY, 

yn+i = 2(f)o{u)yn-yn-i- zb'^Y, v = uh, z = eh'^. 
Elimination of the vector Y gives the recurrence relation 

yn+i - ^(z/^ z) yn + P(z/^ z) yn-i = 0, (5.44) 

where 

S{iy\z) =2(f)o{iy) - zb^ N-^ (e + c), P{u\ z) = 1 - zb^ c, (5.45) 

and 

N = I + {u^ + z)A, e={l,...,lf. (5.46) 
The characteristic equation is 

e-S{u\z)^ + P{u\z)=0. (5.47) 
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Firstly, let us consider dissipative ATSH methods. Working with fl5.47p . we can ask, 
for a given method (i.e., a given u), and a given test frequency A, what restriction must 
be placed on the stepsize h to ensure that the stability condition 

P{u^,z)<l and \S{u^,z)\ < P{u^,z) + 1, (5.48) 

is satisfied. This question can be answered by examining S'(i/^, z) and P(z/^, z) in the u— 
z plane. For ARKN methods such a stability analysis was introduced by Franco (2005). 
The following definition is originally formulated by Coleman & Ixaru (1996) for expo- 
nentially fitted methods for (II. ip . Here, it is adjusted in terms of the methods of 
concern. 

Definition 5 For a dissipative ATSH method with S'(z/^, z) and P(z^^, z) where v = ujh 
and z = eh, and uj and e are given, the primary interval of absolute stability is the 
largest interval (0, /iq) such that Iji5.48\ ) holds for all stepsizes h G (0, /iq)- If, when Hq 
is finite, Iji5.48\ ) holds also for 'y < h < 6, where 7 > /iq then the interval (7, 6) is a 
secondary interval of absolute stability. The region of absolute stability is a region in 
the v — z plane (u > 0), throughout which ^5.48\ ) holds. Any closed curve defined by 



P{v\z) = l or \S{v'',z)\=P{v\z) + l, 
is a stability boundary. 

Likewise, for zero-dissipative ATSH methods the definition of the primary interval of 
periodicity and the region of periodicity is evident. 

In the particular case when the main frequency is exactly known (i.e. 2; = 0) we 
have for both dissipative and zero-dissipative methods that 

5(i/^ 0) = 2 cos(z/) and P(z/2, 0) = 1. 

It follows that the z/-axis is a stability boundary. On this line the periodicity condi- 
tion (12.131) is satisfied except when v = nix ioi positive integer n. 

In the dissipative case, when the frequency is not exactly known the stepsize has 
to be selected carefully. Here we show some sensible points. 

Theorem 3 For dissipative ATSH methods there exist values for uj and e for which 
the primary interval of absolute stability is empty. 
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Proof. Consider the function F defined as 

F{H^) = {I + A)~^ c. 

Assume that F is continuous at = . So we can find an interval (— -Zq, zq) such that 
FiiP' + 2;) has the same sign for all z G (—2:0,-20). It turns out that for such ^-values 
the function P, as given in fl5.45p - fl5.46p . has a different sign at the points (z/, —2;) and 
(z/, 2;). From the absolute stability condition 05.481) it follows that an ATSH method 
which is stable at (z/, — z), is not stable at (z/, z). Thus the z/-axis acts as a stability 
boundary in the sense that it separates stable and unstable regions. This concludes 
the proof. 



6 Phase-lag and dissipation analysis 

For any method corresponding to the characteristic equation f l5.47p . the quantities 

^{v\ z)=H- arccos f , d(z/^ z) = I - y^P{^, (6.49) 

\2^P{u^,z)J 

are called the phase-lag and the amplification error, respectively. As pointed out by 
Franco (2005) for ARKN methods, the analysis of the phase-lag and the dissipation 
becomes more useful if we introduce 

H, z = -^H\ (6.50) 



in fl6.49p . So we arrive to the following definition. 

Definition 6 The phase-lag order is q if 

cl>[v\ z) = c^{u\ e) H'^+' + 0(if''+3), (6.51) 
and the dissipation order is r if 

d{u^, z) = Cd{uj\ e) iy'+^ + 0{H'+^). (6.52) 

Ctf){uj'^,e) and Cd(u;^,e) are called the phase-lag and dissipation functions, respectively. 

In the particular case when the main frequency is exactly known (i.e. 2; = 0) the test 
equation (15. 430 is integrated exactly and so there is no phase-error and no dissipation. 
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We investigate the phase properties when the main frequency is not exactly 
known. Let us define Cj := A^~^ c and Uj := b'^ A^~^ e. Some algebraic manipulation 
gives 

• ATSH method of algebraic order p = 2k: 

oo 

+ ^ 5^(-ly(f/,+C,)i^^^ (6.53) 

j=k+l 

oo 

j=k+l 

• ATSH method of algebraic order p = 2 k — 1: 

j=o ^ ' j=k+l ^ ' ^ 

oo 



oo 

e 



j=k+i 



(6.54) 



When substituting ( 16.53p - (16.54p in (16.49P and then considering the Taylor expansion 
with respect to H it is sufficient to retain the term with the lowest power. After tedious 
but straightforward calculations we have concluded with 

Theorem 4 1. Assume that the algebraic order p of a dissipative TSH method is 
even (odd) and that the phase-lag order is q = p (q = p + 1). Then the cor- 
responding ATSH method has also phase-lag order q. The leading term of the 
phase-lag ^6. 51\) is 

c<t>{^, e) = c^, (6.55) 

where is the phase-lag constant of the classical TSH method. 

2. A dissipative TSH method and the corresponding ATSH method have both the 
same dissipation order. The leading term of the dissipation ^6. 52) is 

d^{uj,e) = ^- — d^, (6.56) 
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where Cd is the dissipation constant of the classical TSH method. 



From (16.551) it follows that the conditions for a ATSH method to have phase-lag order 
g = p + 2 (p: even) or g = p + 3 (p: odd) are exactly the same as those of the 
corresponding classical method. This establishes 

Corollary 5 Assume that the algebraic order p of a TSH method is even (odd) and 
that the phase-lag order is q = p + 2 (q = p + 3). Then the corresponding ATSH method 
has also phase-lag order q. 

In general, Scheifele's adaptation does not conserve the phase-lag order for dissipative 
TSH methods. In contrast, we will show that the phase-lag order is always conserved 
in the zero-dissipative case. Taking into account the order conditions obtained in Sec- 
tion [Hand proceeding as in Section 9 of Coleman (2003) we can reformulate Coleman's 
Theorem 6 for zero-dissipative ATSH methods as follows. 

Theorem 6 For the determination of the the phase-lag order of a zero-dissipative 
ATSH method liS. 28\) we have to compute the scalar quantities Ck = 6'^v4^~^c and 
Uk = A^~^ e for k = 1,2,.... The phase-lag order is q iff Uk = 2 (f)2k{^) for k = 
1, . . . , [^^] and Ck = for A; = 1, . . . , [|] but one of those conditions is not satisfied 
when p is replaced by p -\- 1. 

Corollary 7 A zero-dissipative ATSH method and its classical companion have both 
the same phase-lag order. 

The phase-lag function is also of the form (I6.55p . 

Obviously we have in all cases that c<^(c<j,0) = Cd{uj,0) = 0, 0^(0, e) = and 
q(0, e) = Q. When an acceptable estimate of the dominant frequency is available 
(i.e. e ~ 0) the magnitude of the phase-lag (I6.55P and the amplification error (16.561) is 
then much smaller than those of the corresponding classical method. Furthermore, the 
more accurate the estimate of the dominant frequency, the smaller the phase-lag and 
the amplification error. 
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7 Construction of explicit ATSH methods 



In this section we study the construction of exphcit ATSH methods with algebraic 
orders four and five. Both dissipative and zero-dissipative methods are presented. The 
construction procedure in the classical case was previously considered by Franco (2006a). 



7.1 Methods using two function evaluations per step 

Consider the explicit ATSH method defined by the table of coefficients 



-1 

























asi 


032 







bi 


&2 


&3 



Under the simplifying assumptions (see Coleman (2003)) 

Ae='-^, (7.57) 

the order conditions up to order four are 

fe^e = 202('^), b^c = 0, fe^c^ = 4 04(z/), 6^c^ = 0, b^Ac=0. 

(7.58) 

We have the unique solution 

61 = 63 = 204(1^), 62 = -4 04(z/) + 2 02(z^), C3 = 1, a3i = 0, 032 = 1. 

(7.59) 

When u ^ the method reduces to the explicit Numerov method of Chawla (1984). 
Remark that the values (17.591) are obtained in a different way by Van de Vyver (2007a). 
A stability and phase-lag analysis is also included in that paper. 



7.2 Methods using three function evaluations per step 

Next, we analyze the construction of explicit ATSH methods defined by the table of 
coefficients 
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7.2.1 Dissipative fifth-order methods 



The order conditions up to order five are given by fl7.57l) - (l7.58l) with in addition 

6^c^ = 48 06(i^), b"" (c.Ac) = + b^ A = A . (7.61) 

Solving the equations fl7.57p -( 17.58p and f l7.6ip . the coefficients (I7.60p are determined 
in terms of the arbitrary parameter C3. Two different strategies will be described in 
order to get an optimal method. A first option is to determine C3 so that the error 
constant E^'^^^ (14.420 is as small as possible. The second option is to choose C3 so 
that the method has phase-lag order eight. 

* ATSH method with minimized error constant 

When minimizing the error constant Eq^^^ , we obtain a value for C3 which is very 
close (within a distance < 10^"^) to those of a classical method of Franco (2006a), 
C3 = 63/100. For this reason we adopt Franco's method and we conclude with the 
coefficients 



126651 _ 900249 _ 100 5i ^2 (720000 0^ - 124158 06 04 + 6031 0l) 



2000000' 2000000' 305488243 0^ 

_ ^1 ^2 (-8000000 (pl + 886200 06 04 + 2849 0^) _ 20000 Si S2 S3 06 

~ 1311912701 ' ~ 213841770101 ' 

6 (40000 06 - 1323 04) 04 
' ~ 1635^ ' 

2 (15338 01 - 240000 06 04 - 3969 04 02 + 75600 02 06) 



189 5. 



2 

_ 400000000(12 06- 04)04 , _ 3748322 0^ _ 6 _ 3 ^2 

^ ~ 30807^3 ' ^ ~ 9 Si S2 S3 ' ~ Too' ~ 3704' 

5i = 600 06 - 13 04, ^2 = 400 06 - 2104, 53 = 40000 06 -287704. 

(7.62) 

The region of absolute stability is drawn in Figure [H The expressions for the phase- lag 
and dissipation associated to this method are given by 

^' - 3780051^ ^' + ^) - ^mm^) + 

* ATSH Method with phase-lag order eight 

Following Corollary the condition that imposes phase-lag order eight is the same as 
that for the classical method. In the classical case, phase-lag order eight is achieved 
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Figure 1: u — z plot for ATSH method fl7.62p . 



when C3 = 25/28, see Franco (2006a). Guided by Franco's method, we conclude with 
the coefficients 

1325 _ 35775 _ 28 Si S2 (18816 (pj - 2186 06 04 + 53 0^) 

~ 43904' ~ 43904' ~ 42930| ' 



Si S2 (526848 (pl - 51800 06 04 + 475 0^) _ 1568 Si S2 S3 

> — — " 

2( 

hi 



2025 01 ' 107325 01 ' 

2 (9408 06 - 625 04) 04 



53 52 

2 (1418 01 - 625 04 02 - 18816 06 04 + 8400 02 06) 

255^ 

2458624(12 06 - 04)04 _ 162 01 _ 25 _ Si 

Oa 



1325^3 ' " S1S2S3' ' 28' " 304^ 

Si = 336 06 - 25 04, ^2 = 168 06 - 11 04, ^3 = 9408 06 - 775 04- 



(7.63) 
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The region of absolute stability is drawn in Figure [2l The phase-lag and dissipation 
for this method are 
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Figure 2: p — z plot for ATSH method fl7.63p . 



7.2.2 Zero-dissipative fourth-order method with phase-lag order six 

Here we investigate how we can obtain zero-dissipative methods. Following Theorem [6] 
the method has phase- lag order six when 

b^A^c = Q, A^e = 2(t)e{iy). (7.64) 

We find C3 = 1 which is incompatible with the fifth-order conditions fl7.6ip . and the 
algebraic order of the method should be restricted to four. Solving equations (17.571) . 
(I7.58P and (17.640 we obtain the coefficients in terms of arbitrary parameters C3 and 
C4. The error constant E^"^^^ (14.420 should be as small as possible so we have that 
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C4 = (5c3 — 2)/(5c3 — 5), just like Franco's original case. It is easy to verify that the 
method reaches order five for linear systems of ODEs 



-uj'^y + g{x). (7.65) 

In the classical case the free parameter C3 is chosen so that the resulting method is op- 
timized for the class of linear problems fl7.65p . Here, in order to calculate the error con- 
stant when solving (17.651) . we have to consider the coefficients of the 7th-order elemen- 
tary differentials /(5)(a;)(2/',y',y',y',y'), f'Kv) {f^Kx){y' ,y\y')) ^ndu^ f^\x){y' ,y' ,y') 
The other 7th-order elementary differentials remain zero for (I7.65p . Minimizing this er- 
ror constant we obtain C3 = 13/20. For comparison, in the classical case Franco (2006a) 
obtained C3 = 33/50. The following coefficients are found 

5(7640 06 + 63704) 





= 0, a32 = 


429 
" 800' 


O-Al — 


38200 06 
79233 04' 


^42 ■ 




764000 06 


hi = 


6 04 


h2 = - 


596 04 


1030029 04' 


11 


65 


64 = 


4802 04 
955 ' 


13 
~ 20' 


C4 


5 

~ ~7' 





+ 2, 



31213 04 

128000 04 

03 = , 

^ 27313 ' 



(7.66) 



The region of periodicity is drawn in Figure [3], and the phase-lag is 

' ' 40320 u;2 + e ^ ' 



8 Numerical experiments 

In order to evaluate the effectiveness of the new method derived above we consider 
several model problems. The new method have been compared with other explicit TSH 
codes proposed in the literature. The criterion used in the numerical comparisons is the 
usual test based on computing the maximum global error over the whole integration 
interval. In Figures HHS] we have depicted the efficiency curves for the tested codes. 
These figures show the decimal logarithm of the maximum global error versus the 
computational error measured by the number of function evaluations required by each 
code. The algorithms used in the comparisons have been denoted by 
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V 

Figure 3: u — z plot for ATSH method fl7.66p . 

• CHARA6(8,oo): Zero-dissipative method derived by Chawla & Rao (1987). 

• FRA5(8,5): Classical method derived by Franco (2006a). 

• FTSH5(6,5): Phase-fitted and amplification- fitted method derived by Van de 
Vyver (2006). 

• ATSH5(6,5): ATSH method i^M>- 

• ATSH5(8,5): ATSH method (USSD- 

• ATSH4(6,oo): ATSH method (17:^ . 

Here, A(B,C) means that the method has algebraic order A, phase-lag order B and 
dissipation order C. 

We have used the following five model problems: 

Problem 1. An inhomogeneous equation studied by van der Houwen and Som- 
meijer (1987) 

/ = -100?/ + 99 sin(x), 1/(0) = 1, y'(0) = 11. 
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The exact solution is given by: 

y{x) = cos(10 x) + sin(10 x) + sin(a;). 

It consists of a rapidly and slowly oscillating function; the slowly varying function is due 
to the inhomogeneous term. The equation has been solved in the interval [0, 100] with 
fitted frequency is c<j = 10. The numerical results stated in Fig. Hlhave been computed 
with stepsizes h = 2"^ j = 2, . . . , 6 for CHARA6(8,cx)) and FTSH5(6,5), j = 3, . . . , 7 
for FRA5(8,5) and j = 1, . . . , 5 for ATSH5(6,5), ATSH5(8,5) and ATSH4(6,oo). 

Problem 2. An "almost periodic" orbit problem studied by Stiefel and Bet- 
tis (1969) 

^" = -^ + 0.001 e*"^, ^(0) = 1, ;z'(0) = 0.9995 i 

The equation has been solved in the interval [0, 1000] with fitted frequency uj = 1. The 
exact solution is given by: 

z{x) = (1 - 0.0005 ix) e'"". 

The solution represents a motion of a perturbation of a circular orbit in the complex 
plane. The problem may be solved either as a single equation in complex arithmetic 
or as a pair of uncoupled equations. The numerical results stated in Fig. H] have been 
computed with stepsizes h = 2~\ j = -2,..., 2 for CHARA6(8,oo), ATSH5(6,5) 
and ATSH4(6,cx)), j = -1,...,3 for FTSH5(6,5) and ATSH5(8,5), j = 0,...,4 for 
FRA5(8,5). 

Problem 3. A satellite problem studied by Ferrdndiz et al. (1992) 
We consider the problem of determining the position of an earth satellite. The equa- 
tions of motion have been expressed in focal variables (see Ferrandiz (1988) and 
Ferrandiz et al. (1992)). The coordinates of the basic set of focal variables are three 
components (?/i, ?/2; l/s) of the direction vector of the particle and the inverse u of the 
radial distance. In this formulation the satellite problem can be formulated in four 
decoupled pertubed harmonic oscillators with unit frequency: 

y'l + yi = Qh i = 1, 2, 3, 

u" + u = ^ + Q, 

where /i is the reduced mass, while Qi and Q denote the corresponding perturbation 
terms. We consider the almost periodic equatorial orbit with the zonal harmonic 
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coefficient J2 taken as the perturbation parameter. We have neglected higher order 
terms of J2. The system of equations (18.671) can be written in the form 

y'l + Vi = 0, ^ = 1,2,3, 

J ^^-^^^ 

u" + u = 4 + 124^', 
c c 

where c is the angular momentum and it can be considered as a constant. The solutions 
of the ffist three oscillators are trivial thus we are focused on the last equation. We 
consider the domain of integration [vr, 100] . The initial conditions are given by 

, , u (1 — e) , , , 

For our numerical purpose we consider orbits with eccentricity e = 0.99. In this case: 

_ 100 J2 _ 50 
c2 ~ 20895' c2 ~ 20895000' 

The error has been calculated using a reference solution obtained by means of the 
perturbation techniques developed by Farto et al. (1998). The numerical results stated 
in Fig. [5] have been computed with stepsizes h = {1 — 7r/100)2~-^, j = — 1,...,3 
for CHARA6(8,oo) and FTSH5(6,5), j = 0,...,4 for FRA5(8,5), j = -2, . . . , 2 for 
ATSH5(6,5), ATSH5(8,5) and ATSH4(6,oo). 

Problem 4. A perturbed system studied by Franco (2002) 
As an example of a system we consider 

y'i = -25 - e {yl + yl) + th{x), yi(0) = 1, y'M = 0, 
y'> = _25y2 - e (y? + yl) + e/2(x), ^2(0) = e, y'M = 5, 

where 

/i(x) = 1 + + 2e sin(5x + a;2) + 2 cos(x2) + (25 - 4x2) sin(x2), 
/2(x) = 1 + e2 + 2e sin(5x + x2) - 2 sin(x2) + (25 - 4x2) cos(x2). 

In our test we choose e = 10"^. The system has been solved in the interval [0, 5] with 
= 5. The analytical solution is given by: 

yi{x) = cos(5x) + e sin(x^), ^2(2;) = sin(5x) + e cos(a;^). 

The numerical results stated in Fig. \5\ have been computed with stepsizes h = 2~\ 
j = 1, . . . , 5 for CHARA6(8,oo), j = 2, . . . , 6 for the other codes. 
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9 Conclusions 



Scheifele's G-functions methods are designed in such a way that the exact integration of 
the homogeneous solution of perturbed oscillators fll.2p is automatically included. The 
methods take care with the evaluation of the inhomogeneous part of (11.21) . i.e. g{x,y). 
We have applied Scheifele's approach to TSH methods for an accurate and efficient 
integration of (11. 2p . The resulting methods, called ATSH methods, have coefficients 
dependent onu = u h, where u is a specified angular frequency. Classical TSH methods 
are the limiting forms of ATSH methods as z/ — 0. 

This paper provides a theoretical framework for the derivation of ATSH methods. 
One of our main aims is to develop the order conditions for this new type of methods. It 
is found that ATSH methods share some important properties with the corresponding 
classical TSH methods such as zero-stability, the dissipation order and, under some 
conditions, with the phase-lag order. On the contrary, the stability properties are very 
different from the classical method and they depend on the fitted frequency and the 
stepsize. When the main frequency of the problem is exactly known stability problems 
will never occur, except for a discrete set of exceptional values of the stepsize. When 
the dominant frequency is not exactly known some care is required when selecting the 
stepsize. 

In particular, we have demonstrated the validity of the theory with explicit 
fourth- and fifth-order ATSH methods. The new methods are adaptations of the clas- 
sical TSH methods of Franco (2006a). In most cases, the dissipative ATSH method 
(I7.62P with minimized error constant outperforms all the other methods considered. In 
contrast with the results of the phase-fitted and amplification-fitted methods of Van de 
Vyver (2007b), it turns out that the accuracy of ATSH methods is mostly determined 
by its usual local truncation error rather than by its phase-lag. 

Our task is restricted to scalar equations or systems involving only one frequency. 
When solving systems with more than one frequency, or more general, systems of the 
form 

y" = Ky + g{x,y), (9.69) 

the resulting methods have coefficients which are functions of the matrix h"^ K. So their 
evaluation is not direct. To overcome this difficulty, together with some other troubles. 
Franco (2006b) has modified ARKN methods for oscillatory systems of the form (19.690 . 
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The extension of Franco's approach to the ATSH methods considered here might be 
an interesting suggestion for some future work. 
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Figure 4: Efficiency curves of tfie metfiods for Problems 1-2. 
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